#ifndef __NAVIERSTOKES_H
#define __NAVIERSTOKES_H 1

#include <feel/options.hpp>
#include <feel/feelcore/application.hpp>

#include <feel/feelalg/backend.hpp>
#include <feel/feeldiscr/functionspace.hpp>
#include <feel/feelfilters/gmsh.hpp>
#include <feel/feelfilters/exporter.hpp>

#include <feel/feelmesh/elements.hpp>
#include <feel/feelvf/vf.hpp>
#include <feel/feeldiscr/bdf2.hpp>
#include <feel/feeldiscr/ale.hpp>
#include <feel/feeldiscr/operatorinterpolation.hpp>

namespace Feel
{
    namespace FSI
    {

    template<int Dim,
             int Order,
             template<uint16_type,uint16_type,uint16_type> class Entity>
    class NavierStokes : public boost::enable_shared_from_this< NavierStokes<Dim, Order, Entity> >
    {
    public:

        typedef NavierStokes<Dim, Order, Entity> self_type;
        typedef boost::shared_ptr<self_type> self_ptrtype;

        //___________________________________________________________________________________//

        typedef double value_type;

        typedef Application application_type;
        typedef boost::shared_ptr<application_type> application_ptrtype;

        //___________________________________________________________________________________//

        typedef Backend<value_type> backend_type;
        typedef boost::shared_ptr<backend_type> backend_ptrtype;

        //___________________________________________________________________________________//

        /*matrix*/
        typedef typename backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
        typedef typename backend_type::vector_ptrtype vector_ptrtype;

        //___________________________________________________________________________________//

        /*mesh*/
        typedef Entity<Dim,1,Dim> convex_type;
        typedef Mesh<convex_type> mesh_type;
        typedef boost::shared_ptr<mesh_type> mesh_ptrtype;

        //___________________________________________________________________________________//

        /*basis fluid*/
        typedef Lagrange<Order, Vectorial,Continuous,PointSetFekete> basis_fluid_u_type;
        typedef Lagrange<Order-1, Scalar,Continuous,PointSetFekete> basis_fluid_p_type;
        typedef Lagrange<0, Scalar> basis_l_type;

        typedef bases<basis_fluid_u_type,basis_fluid_p_type> basis_fluid_type;

        /*space fluid*/
        typedef FunctionSpace<mesh_type, basis_fluid_type> space_fluid_type;
        typedef boost::shared_ptr<space_fluid_type> space_fluid_ptrtype;

        BOOST_MPL_ASSERT( ( boost::is_same<typename space_fluid_type::bases_list, basis_fluid_type> ) );
        BOOST_MPL_ASSERT( ( boost::is_same<typename mpl::at<typename space_fluid_type::bases_list,mpl::int_<0> >::type, basis_fluid_u_type> ) );
        BOOST_MPL_ASSERT( ( boost::is_same<typename mpl::at<typename space_fluid_type::bases_list,mpl::int_<1> >::type, basis_fluid_p_type> ) );

        /* functions */
        typedef typename space_fluid_type::element_type element_fluid_type;
        typedef boost::shared_ptr<element_fluid_type> element_fluid_ptrtype;

        typedef typename element_fluid_type::template sub_element<0>::type element_fluid_0_type;
        typedef typename element_fluid_type::template sub_element<1>::type element_fluid_1_type;

        //___________________________________________________________________________________//

        /* time */
        typedef Bdf<space_fluid_type>  bdf_type;
        typedef boost::shared_ptr<bdf_type> bdf_ptrtype;

        //___________________________________________________________________________________//

        /* export */
        typedef Exporter<mesh_type> export_type;

        //___________________________________________________________________________________//
        typedef typename space_fluid_type::template sub_functionspace<1>::type::element_type space_fluid_pressure_type;
        typedef typename space_fluid_type::template sub_functionspace<1>::type space_fluid_pressure_ptrtype;
        typedef typename space_fluid_pressure_type::element_type element_fluid_pressure_type;
        typedef boost::shared_ptr<element_fluid_pressure_type> element_fluid_pressure_ptrtype;

        /* ALE */ /*verifier l'utilite de toutes les def!"*/
        typedef typename space_fluid_type::template sub_functionspace<0>::type::element_type space_fluid_velocity_type;
        typedef typename space_fluid_type::template sub_functionspace<0>::type space_fluid_velocity_ptrtype;
        typedef typename space_fluid_velocity_type::element_type element_fluid_velocity_type;
        typedef boost::shared_ptr<element_fluid_velocity_type> element_fluid_velocity_ptrtype;

        typedef Bdf< space_fluid_velocity_type > bdf_ale_velocity_type;
        typedef boost::shared_ptr<bdf_ale_velocity_type> bdf_ale_velocity_ptrtype;

        typedef ALE< convex_type > ale_map_type;
        typedef boost::shared_ptr< ale_map_type > ale_map_ptrtype;

        typedef typename ale_map_type::pN_functionspace_type ale_map_functionspace_type;
        typedef typename ale_map_type::pN_element_type ale_map_element_type;
        typedef boost::shared_ptr<ale_map_element_type> ale_map_element_ptrtype;

        typedef Bdf< ale_map_functionspace_type > bdf_ale_map_type;
        typedef boost::shared_ptr<bdf_ale_map_type> bdf_ale_map_ptrtype;

        // Interface
        typedef Mesh< Simplex<1, 1,2> > mesh_structure_type;
        typedef boost::shared_ptr<mesh_structure_type> mesh_structure_ptrtype;

        typedef Lagrange< Order, Scalar/*, Continuous, double, PointSetWarpBlend*/ > basis_structure_type;
        typedef FunctionSpace<mesh_structure_type, basis_structure_type > structure_space_type;
        typedef boost::shared_ptr<structure_space_type> structure_space_ptrtype;
        typedef typename structure_space_type::element_type element_structure_type;
        typedef boost::shared_ptr<element_structure_type> element_structure_ptrtype;

        //___________________________________________________________________________________//

        typedef element_fluid_0_type element_dirichlet_type;

        //___________________________________________________________________________________//

        mesh_ptrtype createFluidMesh( double meshSize, std::string name );

        mesh_structure_ptrtype createStructMesh( double meshSize, double __ypos, std::string name );

        //___________________________________________________________________________________//

        NavierStokes( application_ptrtype __appli );

        NavierStokes( NavierStokes const & M ) { std::cout <<"\nWARNING constructeur copi\n"; }

        //___________________________________________________________________________________//

        boost::shared_ptr<element_dirichlet_type> getDirichletElement()
        {
            element_fluid_type bid( M_Xh, "bid" );
            boost::shared_ptr<element_fluid_0_type> bid_u(new element_fluid_0_type(bid.template element<0>()));
            return bid_u;
        }

        mesh_ptrtype mesh() { return mesh_fluid; }

        space_fluid_ptrtype functionSpace() { return M_Xh; }

        backend_ptrtype backend() { M_backend; }

        element_fluid_ptrtype getSolution() { return U_fluid; }

        element_fluid_velocity_type & meshVelocity() { return *M_mesh_velocity; }

        element_fluid_ptrtype getNormalStress();

        bool weakCL() { return M_weakCL; }

        void changeRepository(std::string pathLoc="");

        //___________________________________________________________________________________//

        void runPtFixe();

        void solveFluid(std::string __PDEtype);

        //___________________________________________________________________________________//

        void init();
        void updateJacobian( const vector_ptrtype& X, sparse_matrix_ptrtype& J , vector_ptrtype& R);
        void updateResidual( const vector_ptrtype& X, vector_ptrtype& R );


        void updateLinearPDE(const vector_ptrtype& X,sparse_matrix_ptrtype& A , vector_ptrtype& F);


        void updateStokes(/*element_fluid_type & U*/const vector_ptrtype& X, sparse_matrix_ptrtype& A , vector_ptrtype& F);

        void updatePtFixe( element_fluid_type & U_old, sparse_matrix_ptrtype& A , vector_ptrtype& F);

        //___________________________________________________________________________________//


        application_ptrtype application() { return M_lifeApplication; }


        element_structure_ptrtype getDeplTop() { return M_depl_top; }
        element_structure_ptrtype getDeplBottom() { return M_depl_bottom; }

        //___________________________________________________________________________________//

        void updateTime(double __t);

        void updateDeplacement(element_structure_ptrtype __depl_top,
                               element_structure_ptrtype __depl_bottom );


        ////
#if 0
        typedef Feel::Simplex<2,1,2> convex_structStVen_type;
        typedef Mesh<convex_structStVen_type> mesh_structStVen_type;
        typedef boost::shared_ptr<mesh_structStVen_type> mesh_structStVen_ptrtype;

        /*basis*/
        typedef bases<Lagrange<2, Vectorial> > basis_structStVen_type;

        /*space*/
        typedef FunctionSpace<mesh_structStVen_type, basis_structStVen_type,Continuous, double> space_structStVen_type;
        typedef boost::shared_ptr<space_structStVen_type> space_structStVen_ptrtype;

        /* functions */
        typedef typename space_structStVen_type::element_type element_structStVen_type;
        typedef boost::shared_ptr<element_structStVen_type> element_structStVen_ptrtype;

        //typedef typename Feel::StVenantKirchhoff<2,2,Feel::Simplex>::element_struct_type element_structStVen_type;
        //typedef typename Feel::StVenantKirchhoff<2,2,Feel::Simplex>::element_struct_ptrtype element_structStVen_ptrtype;
#endif

        template <typename element_mecasol_ptrtype>
        void updateStructDeplacement(element_mecasol_ptrtype structSol);

        void updateALEmap();

        void updateBdf();

        //___________________________________________________________________________________//

        void exportResults() { exportResultsInstat(mesh_fluid,*U_fluid,*M_mesh_velocity,*M_aux, M_cpt_export); }


private:

        void solve( sparse_matrix_ptrtype const& D, element_fluid_type& u, vector_ptrtype const& F, bool is_sym );

        double diffL2( element_fluid_type & u1, element_fluid_type & u2 );

        void exportResults( element_fluid_type& u );

        void exportResultsInstat( element_fluid_type& u,
                                  element_fluid_velocity_type & mesh_velocity,
                                  ale_map_element_type & aux,
                                  double time );

        void exportResultsInstat( mesh_ptrtype __mesh,
                                  element_fluid_type& u,
                                  element_fluid_velocity_type & mesh_velocity,
                                  ale_map_element_type & aux,
                                  double time );


    private:

        application_ptrtype M_lifeApplication;

        backend_ptrtype M_backend;

        double meshSize;

        mesh_ptrtype mesh_fluid, mesh_reference_fluid;
        MeshMover<mesh_type> M_mesh_mover;

        space_fluid_ptrtype M_Xh,M_Xh_ref;
        element_fluid_ptrtype U_fluid,M_Ufluid_ref;

        element_fluid_velocity_ptrtype M_identity_ale;
        element_fluid_velocity_ptrtype M_mesh_velocity;

        bdf_ptrtype M_bdf_fluid;


        mesh_structure_ptrtype mesh_struct_top;
        mesh_structure_ptrtype mesh_struct_bottom;
        structure_space_ptrtype Sh_top;
        structure_space_ptrtype Sh_bottom;
        element_structure_ptrtype M_depl_top;
        element_structure_ptrtype M_depl_bottom;

        element_structure_ptrtype M_bc_reference_top;
        element_structure_ptrtype M_bc_reference_bottom;

        std::vector<element_structure_type> M_referencePolyBoundarySet;

        boost::shared_ptr< ALE< convex_type > > M_aleFactory;
        ale_map_element_ptrtype M_aux;

        bdf_ale_map_ptrtype M_bdf_ale_map;
        bdf_ale_velocity_ptrtype M_bdf_ale_identity;

        double M_time;
        uint M_cpt_export;
        double rho;
        double M_mu;
        double penalbc;
        double Vitesse_Entree;
        double Pression_Sortie;

        boost::shared_ptr<export_type> M_exporter;
        boost::shared_ptr<export_type> M_exporter_ref;

        bool M_weakCL;
        bool M_isStat;
        bool M_isMoveDomain;
        bool M_isFixeBoundary;

        bool M_LookComment;



    }; // Stokes







template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
template <typename element_mecasol_ptrtype>
void
NavierStokes<Dim, Order, Entity>::updateStructDeplacement( element_mecasol_ptrtype structSol )
{

    if(M_LookComment) std::cout<<"[NavierStokes] : updateStructDeplacement start\n";

    this->application()->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
                                           % this->application()->about().appName()
                                           % convex_type::name()
                                           % Order
                                           % this->meshSize
                                           );

    typedef typename element_mecasol_ptrtype::element_type element_mecasol_type;
    typedef typename element_mecasol_type::functionspace_type space_mecasol_type;

    element_mecasol_type uStruct(structSol->functionSpace(),"jujuj");
    uStruct = vf::project( uStruct.functionSpace(),
                           elements( uStruct.functionSpace()->mesh()),
                           vf::trans(vf::idv(*structSol))*vf::oneY()*vf::one() );
                           //vf::trans(vf::idv(*structSol))*vf::N()*vf::one() );

    //interpolate on each moving boundary
    //OperatorInterpolation<space_structStVen_type, structure_space_type > opIbottom( structSol->functionSpace(), Sh_bottom , M_backend );
    OperatorInterpolation<space_mecasol_type, structure_space_type > opIbottom( structSol->functionSpace(), Sh_bottom , M_backend );
    opIbottom.apply(uStruct,*M_depl_bottom);

    //OperatorInterpolation<space_structStVen_type, structure_space_type > opItop( structSol->functionSpace(), Sh_top , M_backend );
    OperatorInterpolation<space_mecasol_type, structure_space_type > opItop( structSol->functionSpace(), Sh_top , M_backend );
    opItop.apply(uStruct,*M_depl_top);

    //std::cout << "\n norm l2 deplacement bottom : " << M_depl_bottom->l2Norm()<<"\n";
    //std::cout << "\n norm l2 deplacement top : " << M_depl_top->l2Norm()<<"\n";

    if (M_isMoveDomain)
        {
            this->updateALEmap();
        }

    if (M_LookComment) std::cout<<"[NavierStokes] : updateStructDeplacement finish\n";


}



    } // NavierStokes

} // Feel






#endif /* __NSSTOKES_H */

